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1 Introduction 

The problem of concern here is the analysis of partially polarized light re- 
flected by plant canopies. The significance of such signals and certain under- 
lying mathematical models are presented in [3], [4], [5], for example. 

We start with data measurements obtained at specified grid points with 
instruments with a particular field of view (FOV). The goal is to generate 
equivalent data values which would correspond to a more desirable array 
of grid points, such as uniformly spaced ones, and possibly with a different 
FOV. Then, polarization parameters, specifically the polarized reflectance 
Rqu and angle of the polarization plane x, are computed. These quantities 
are related to the quantities Rq and Ru which are linearly related to raw 
sensor measurements as: 


Rque j2x — Rq + jRu ( 1 ) 

The basic approach is to model the measurement process as linear, based 
on underlying “detail”, obtaining an estimate of the detail and recomput- 
ing “simulated” measurements on the new grid with the new FOV. Actu- 
ally, these two stages are combined into one so the dense “detail” is not 
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actually computed. Because of imperfections in the measurement process, 
least-square methods are used. The criterion for matching is a quadratic 
one directly on the data values themselves, or, alternatively, a quadratic one 
based on a nonlinear transformation of the measured data. 

The solution formula in the linear case involves pseudo-inverses of matri- 
ces, and hence is computed using singular- value decomposition (SVD) meth- 
ods. 

The primary physical processes underlying the measurement process are: 

• the blurring effect of a finite FOV; 

• the variability of the area of the ground “footprint” with view angle; 

• leaf size and orientation, ground cover and heliotropy; in particular, 
the solar position is not constant for all data measurements, and he- 
liotropism introduces a “bias” in the probabilistic model of the leaf 
orientations; 

• the different sensors were not collocated, so measurements taken at 
the same time at the same view angle result in partially overlapping 
“footprints” ; 

• the data fusion effect; that is, whether the data from different sen- 
sors, at different times, or at different wavelengths should be treated 
separately or combined, for example averaged; 

• the desired polarization information, such as Rqu, is computed as a two 
step process: the first is a linear transformation of information from 
three different sensors (yielding Rq , Ru) and the second is a nonlinear 

one (Rqu — yj R q + Ru > X — tan -1 ( Ru/Rq ) /2). Whether the linearly 
dependent parameters are recovered, then the nonlinearity is applied, 
or if the nonlinearly dependent parameters are computed directly does 
effect the ultimate outcome. 

Whatever the details in the recovery method, the results must then be 
compared against those predicted by theoretical models. These theoretical 
models ultimately depend on models of probability density functions for var- 
ious parameters, such as leaf orientation and ground cover, and their relation 
to source and view angle. 

In the linear approach, these issues are addressed in the following fashion. 
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the measurement process is a linear one which is relatively localized; 
hence it can be represented as a matrix acting on a data vector. 


• the FOV concept yields a sparse matrix; if a single matrix is applied 
to the entire canopy area, the structure is not simple (for example, 
not banded) because of the way data is organized in the vector, corre- 
sponding to irregular 2-D sampling; alternatively, for every view angle 
for which a computation is be performed, a full, lower dimensional ma- 
trix can be used to model the local measurement process and pseudo- 
inversion methods applied locally; this latter approach can be taken as 
a spatially varying by linear filtering process which must be inverted; 

• data fusion can be handled by either treating data sources separately, 
or combining them into a joint vector; 

• the linear transformation step in computing 1 Rq, Rjj can be absorbed in 
the measurement matrix; however, this requires dealing with the three 
sensors together, and hence can cause inaccuracies because the sensors 
are not collocated. 

In these notes, a linear least-squares method is developed, followed by one 
based on a nonlinear cost function. The methods were implemented using 
MATLAB for Windows on a Pentium PC and applied to data collected from 
a Sunflower field in 1991 [2]. The approach here is an outgrowth of the 
one developed in [1]. Several graphical results, shown in the Appendix, are 
discussed. 


2 Linear Least-Squares Method 

2.1 General Notation 

In what follows, let A f [B) and 7Z (B) denote the nullspace and range space 
of a matrix £, and B& the pseudo- inverse of B. Every matrix B is a one- 
to-one onto linear map from J\f L ( B ) to 7 Z(B), and is the inverse map. 
In particular, Bx € 1Z ( B ) and B#y £ A/* 1 (B) for all x,y. The matrix B#B 
is the orthogonal projection onto A f x (B), and I — B#B is the orthogonal 
projection onto J\f (B). Also, BB# is the orthogonal projection onto 7Z(B). 
Let || -|| denote the L 2 -norm of a vector, and the spectral norm of a matrix. 
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Let B* denote the adjoint of B : the transpose for a real matrix, the Hermitian 
(conjugate) transpose for a complex matrix. 

2.2 The Measurement Process 

We first model the measurement process. By default, all data vectors are 
column vectors and linear operators are represented by matrices acting on 
the left. Let c denote the measured data, obtained at specified grid points 
with instruments with a particular field of view (FOV). Thus, let x represent 
the “underlying detail,” and the smoothing and sampling operations by the 
matrix A, resulting in the measurement process: 

c = Ax An- (2) 

where n is the measurement noise. In general, A is rectangular and has fewer 
rows than columns. Let C be the space of all measurement vectors c, and X 
the space of all detail vectors x. 

Our goal is to compute a data vector, d, corresponding to simulated mea- 
surements onto another set (say, a uniform pattern) of grid points and pos- 
sibly at another FOV. The simulated measurement process is characterized 
by the matrix B. If x were known then: 


d = Bx 


(3) 


Let V denote the space of all simulated measurements d. 

2.3 The Criterion 

We first set up a quadratic criterion which yields a linear solution. Given c, 
the minimum-norm least-square solution for x minimizes || c — Ax ||; since A 
does not have full column rank, in general, there are infinitely many mini- 
mizers, so the specific solution is one which, among all minimizers, also has 
minimum ||aj||. Denoting this solution by x Q , we have: 

x Q = A^c (4) 

and more generally all minimizers of ||c — Ax || are given by: 

x = A#c + x n (5) 
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where x n E J\f (A): 


(6) 


Ax n = 0 

Let Q be the orthogonal projection onto J\f (A): 

Q = I- A* A (7) 

Now, Q # = Q* = Q. Moreover, x n E A f (A) iff it satisfies: 

Xn Qx n (8) 

The choice of x n is arbitrary given only the information in c, and hence 
its effect on computed values constitutes artifacts. 

Simulated measurements d based on such x have the form: 

d = Bx — BA*c + Bx n — BA*c + BQx' (9) 

where x' is arbitrary; that is, given c, then any vector of the form: 

d — BA^c + d Q (10) 

where d a E 1Z ( BQ ), that is d Q which is a linear combination of the columns 
of BQ : corresponds to minimizing ||c — Ax ||. The choice of x n = 0, that is 
x = x Q = A*c achieves minimum ||rc|| but not minimum \\d\\. 

Smoothness (regularity) considerations suggest we may attempt to find 
a solution which minimizes ||d||, not the “underlying detail,” which is never 
directly computed. We can take this one step further. Assume a-priori 
information in the form of a “nominal” cfi is available, so we want to minimize 
\\d-d 1 \\. 

Thus, the problem is as follows: given c and di (which is 0 if there is no 
a-priori information), pick d corresponding to a vector x (via d — Bx) which 
minimizes ||c — Ax\\ and, among all such, \\d — di|| should be minimized. 

Let us denote 7 Z (BQ) as the artifact space. It consists of all vectors in 
the simulated measurement space which arise from measurements performed 
on “underlying detail” which is in the nullspace of the actual measurements. 
Given an actual measurement c, there is no (direct) information on this 
null component, and hence the corresponding component of the simulated 
measurement d can be considered an artifact. 

The orthogonal projector onto the artifact space is P given by: 

P = (BQ) (BQ)* = UU* = jr u k u* k (11) 

k= 1 
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where U = y U\ U 2 ••• u r j is the Lxr matrix whose orthonormal columns, 
the right singular vectors of BQ, span the artifact space. Note that r = 
rank ( BQ ) is the dimension of this artifact space. 

Theorem 1 Given measurements c and a-priori information d\, the simu- 
lated data d which corresponds to minimum \\d — di\\ subject to the constraint 
that d — Bx for some x minimizing ||c — Ax || is given by: 

d=(l-p) BA*c + Pd x (12) 

That is, the artifact component is computed directly from d\, and the compo- 
nent orthogonal to the artifact component is computed directly from c. 

Proof. We have x — AAc + Qx a , with x Q arbitrary, the complete list of x 


which minimize \\c — Ax ||. Then: 


d — Bx = BA^c + BQx a 

(13) 

Then: 


d — di = (yBA^c — di ) + BQx 0 

(14) 

A choice of x a which minimizes: 


d\ — BA^c — BQx 0 

(15) 

is: 

= (BQ)* (d, - BA*c ) 

(16) 

and corresponding d is: 


d = BA^c + Pdi — PBA^c 

(17) 


The vector d in (12) corresponds to detail x which is a least-squares match 
to the measurements c, while minimizing error from a-priori information d\. 
But in the case of no a-priori information, we zero out the artifact component 
via: 

d = (. I-P)BA*c (18) 

= (i-b(i-a*a)(b(i-a*a))*^ba*c 

This involves two pseudo-inverse and five matrix multiplication operations; 
matrix addition and matrix- vector operations are of significantly lower com- 
putational complexity. 

Localizing this reduces complexity to solve for a particular point or region, 
but must be repeated as the matrices need to be reformulated because of the 
irregularity of the original measurements. 
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3 Nonlinear Cost Function 


The underlying cost functions in the previous sections are L 2 - norms, specif- 
ically \\d — di\\. As an alternative, we may try to match results computed 
from d , rather than d itself. For example, suppose we try to match Rqu to 
some a-priori knowledge or probabilistic model. This is a nonlinear function 
of a vector of low-dimension, typically three or less (here, Rqu is a nonlinear 
function of Rq and Ru)- 
Denote this as: 

P (£) (19) 

where: 

£ ~ (fu?2>' ’ * »£at) (20) 

where 1 < N < 3, typically. Hence, d and c are L x N and M x N matrices, 
respectively, with each row (length N) consisting of one vector for which p 
can be computed. Let: 

c = [Cl c 2 • • • Cm J (21) 

1 T 

d = di d 2 • • * dL 

where each c* and di is an A^-dimensional column vector. We set up the cost 
function: 

J (d, di) = ^ [p (di) - pi] 2 (22) 

i - 1 

where p { is p applied to the i th row of d \ . We then find d of the form: 

d = BA^c BQx 0 (23) 

which minimizes (22). Now, we can solve this constrained minimization 
problem using the method of Lagrange multipliers. However, we take an al- 
ternate route here. Instead, we obtain an explicit expression for d satisfying 
the constraints in terms of a minimal set of arbitrary parameters, and sub- 
stitute it into the formula (22) for the cost function J. The unconstrained 
minimization problem can then be solved directly. 

We can rewrite the constraint as: 

(I — P) (cl — BA*c ) = 0 (24) 
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or, in other words, d — BA^c must be in the range space of BQ, that is 
an artifact component. Let {si}[ =1 be orthonormal vectors which span the 


u x u 2 


U r 


artifact space, so that P = Ylk=i u k u t = jS'jS'*, where U = 
is the Lxr matrix whose orthonormal columns span the artifact space. Note 
that r = rank ( BQ ) is the dimension of this artifact space. Hence: 


d = d 0 + u klk (25) 

k = 1 

where: 

d 0 = BA*c (26) 

is L x N is the simulated data computed without regard for artifacts, and 
each 7 fc is 1 x N. The r x N elements in the matrix T given by: 


r 


T T T 

7i 72 ’"Tm 


T 


are free parameters to be chosen so as to minimize J: 


J(?) = E 

2=1 


P U iklkj ~ Pi 


(27) 


(28) 


where is the i th component of the vector d oi is the i th row of the matrix 
d Q and p { is the nominal or a-priori value of p (d) to be matched. 

For 1 < j < N, let 

(I) =dp{u)/dij (29) 

Then the partial derivative of J (T) with respect to the mj th element of 
F is: 


dJ 

dlmj 


2£ 



T u iklk 
k = 1 / 



Pi 


d r 


+ U iklk 
k = 1 


Ur 


(30) 


and the solution to the minimization problem is achieved by setting this 
expression to 0 for 1 < m < iV, 1 < j < r. 

We can simplify the notation at this point somewhat. We can define the 
N x L matrix G via: 


9ji 



+ U ikl k 
k = 1 


Pi 


+ Y; U iklk 

k = 1 


(31) 
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where 1 < j < N, 1 < % < L. Then (30) reduces to: 

{dJ/drf = 2GU (32) 

where (dJ/OTy is the N xr matrix whose jm th element is dJ/d"f m j. Hence, 
the minimization problem is: 

GU = 0 (33) 

This is a system of N • r equations (since GU is an N x r matrix) in N • r 
unknowns, the 7 -’s. 

Some comments at this point is warranted. In general, L may be large, 
corresponding to the number of local or data points being considered. How- 
ever, N is generally small, maybe 2, and r, the rank of the artifact space, is 
selectable. It can be determined either by numerical studies of P, or bounded 
by an arbitrary bound on the computational procedure to be followed. It is 
generally easier to minimize a single nonlinear function of several variables 
rather than solving simultaneous nonlinear equations. We can do this by 
converting the system (33) into a single function: 

N r 

J' (r) = EE ( GU )lk (34) 

m= 1 fc=l 

which is a nonnegative definite form, which achieves its minimum, 0, exactly 
when (33) is satisfied. We can rewrite (34) as: 

J' (T) = trace {GU ( GU )*) = trace (1 GUU*G *) (35) 

and, since P = UU*, we get the final form: 

J'{T) = trace (GPG*) (36) 

We have J' as a nonnegative definite nonlinear quadratic form, and its min- 
imization yields the desired solution. 

The key step here is the choice of r, which is tantamount to a subspace 
selection. 

Example 2 Let the underlying parameters be Rq and Ru, and Rq V = 
Then N = 2, p (£ lt Q = \, and: 

MCi,e 2 ) = Up(Zu&) (37) 

A 6 2 (Cl? £2) = C2/p(Cl?^2) 
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Example 3 If instead we work with the angle of polarization x — tan 1 ( Ru/Rq ) /2. 
Then p(f 1? £ 2 ) = tan -1 (f 2 /fi) / 2 and: 

= — £2/2 (£1 + £2) (38) 

h2 (£n £2) ~ Cl/ 2 (£1 + £2) 

4 Application to Experimental Data and Re- 
sults 

The linear least-squares method and the method based on the nonlinear cost 
function were both applied to data collected from a sunflower field in [[2]]. 

The experiment involved two types of instruments: Barnes, which had a FOV 
of 1°, and Cimel, which had a FOV of 1° or 12°, depending on the day upon 
which measurements were performed. Results for data collected on 07/25/91, 
shown in the Appendix, are representative; on this date, the Cimel FOV was 
12°. 

The measurements were taken from a 13.5m tower in the center of a 
sunflower field approximately 100m x 100m in sunlight. The polarized com- 
ponent of the measured reflectance is modelled [3], [4], [5] as caused by sunlight 
reflected specularly from the leaf surfaces, and hence Fresnel equations for 
specular reflection apply. In particular, we are interested in measuring the 
total polarized reflectance Rqu and the angle of the plane of polarization, y. 

The measurements were taken over multiple zenith angles for several az- 
imuths, approximately ranging from 0° to 360° in increments of 45°. The 
angles of measurement were very irregular. In order to apply the methods 
discussed here, the measurement process had to be modelled by certain op- 
erators, specifically A and B. 

An operator is developed for a fixed ray with a specified view azimuth and 
a zenith angle ranging from 0 to 90°; actually, the zenith angle was restricted 
to a subrange since very small and very large zenith angles caused spots on 
the ground (due to the finite FOV) which either were directly below the tower 
or beyond the edge of the field. Each raw measurement c, underlying x and 
simulated d is obtained from light collected within a finite FOV. The locations 
of the spots at which c is measured is determined by the experimental data 
itself, and is generally irregular; multiple measurements from the same spot 
also can occur. The locations of the underlying detail x were on a dense grid 
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within the vicinity of this strip, and of d on a regular range of spots for the 
fixed azimuth. The azimuth of the x spots varied to encompass the width of 
the spot measurements of c and d. The underlying detail x was assumed to 
arise from a FOV of 1°. 

Given a position of x, it corresponds to a spot whose center, length and 
width depend on the azimuth and zenith values; this spot partially overlaps 
a corresponding spot of c or d, respectively. Hence, each row of A and B, 
respectively, implements a weighted sum of x values to obtain final c and d 
values, respectively, with the weight proportional to the area of the x spot 
which lies within the c (d) spot. Hence, the averaging effect of the FOV 
for data on one grid set is converted to a possibly different FOV on another 
grid set by weighting linearly according to the area overlap of the spots. 
If two points have a large area overlap, the information is expected to be 
more highly correlated than between two points with a smaller area overlap. 
Distance between points is accounted for indirectly, in that spots which are 
further apart have a smaller common area. However, the area overlap method 
is reasonable considering the physical situation, a multitude of plant leaves 
over a ground cover. 

The results, shown in the Appendix, show that the raw measurements 
of R qu and \ are irregular, both for Barnes instruments with FOV 1° and 
Cimel instruments with FOV 12°. The linear least-squares method provides 
generally good matching results, which are significantly smoothed. For the 
case of the Barnes instruments, the Rqu values are somewhat smaller. This 
may reflect averaging over complex values Rq + j Ru , then taking magni- 
tude; phase cancellation due to noise causes diminished Rqu- The nonlinear 
method yields Rqu results similar to the linear method. The key to the non- 
linear method is the a-priori information p; further studies should use the 
mathematical models developed in, for example, [3], [4], [5], and perhaps quan- 
tities other that p = Rqu for better results. The nonlinear method is still 
promising precisely because it more easily accommodates such an approach. 

The values of y obtained for both Barnes and Cimel were very good, 
surprisingly good considering the high sensitivity of phase information in 
many situations. The nonlinear approach produced very small y values, 
suggesting recovered Ru values are too small relative to Rq. This may also be 
caused by the matching to Rqu in the nonlinear method; perhaps combining 
matching to Rqu and x would yield better results. 

In summary, the methods developed here are very promising. Further 
investigations and development of specific implementation in MATLAB, and 
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particularly a more extensive utilization of established mathematical models 
of the underlying physical process, for example through the nonlinear cost 
function approach, is indicated. 
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5 Appendix: Graphs of Results 

Key to linetypes: 

• dots: raw measurements. 

• solid line: linear least-square solution; the a-priori information di is 
taken to be 0. 

• dashed line: solution based on the nonlinear cost function Rqu] the 
a-priori information p is taken to be the mean Rqu corresponding to 
the linear least-square solution do for the specified view azimuth. 
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All graphs are for measurements taken on 07/25/91. 

Barnes instruments had raw FOV of 1° and Cimel instruments had raw 
FOV of 12°, with measurements taken over irregular azimuth and zenith view 
angles. The processed data, by either method, used an effective (‘simulated’) 
FOV of 12° at regularly spaced zenith angles for fixed azimuth angles taken 
from 0° to 315° in increments of 45°. Not all results are shown. 
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